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The elliptical instability is a generic instability which takes place in any rotating flow whose 
streamlines are elliptically deformed. Up to now, it has been widely studied in the case of a constant, 
non-zero differential rotation between the fluid and the elliptical distortion with applications in 
turbulence, aeronautics, planetology and astrophysics. In this letter, we extend previous analytical 
studies and report the first numerical and experimental evidence that elliptical instability can also 
be driven by libration, i.e. periodic oscillations of the differential rotation between the fluid and 
the elliptical distortion, with a zero mean value. Our results suggest that intermittent, space-filling 
turbulence due to this instability can exist in the liquid cores and sub-surface oceans of so-called 
synchronized planets and moons. 

PACS numbers: 47.32.Ef, 95.30.Lz, 47.20.Cq 

The longitudinal libration of a so-called synchronized planet or moon, i.e. the oscillation of its axial rotation rate 
whose mean value is otherwise equal to the orbital rotation rate, arises through its gravitational coupling with its 
closest neighbors [3 [5J. In the body interior, the interaction of a fluid layer (e.g. an iron rich liquid core or a subsurface 
ocean) with the surrounding librating solid shell resulting from viscous, topographic, gravitational or electromagnetic 
coupling, leads to energy dissipation and angular momentum transfer that need to be accounted for in the planet 
thermal history and orbital dynamics, and possibly in its magnetic state [3J. A number of studies has been devoted 
to libration-driven flows in axisymmetric containers in order to investigate the role of the viscous coupling. It has 
been shown that longitudinal libration in an axisymmetric container can drive inertial waves in the bulk of the fluid 
as well as boundary layer centrifugal instabilities in the form of Taylor- Gortler rolls [4j-[9]. In addition, laboratory and 
numerical studies [HHHj have corroborated the analytically predicted generation of a mainly retrograde axisymmetric 
and stationary zonal flow in the bulk, based upon non-linear interactions within the Ekman boundary layers [5J- 
[11] [T3l El ■ Although practical to isolate the effect of viscous coupling, the spherical approximation of the core-mantle 
or ice shell-subsurface ocean boundaries, herein generically called the CMB, is not fully accurate from a planetary 
point of view and very restrictive from a fluid dynamics standpoint. Indeed, due to the rotation of the planet, the 
gravitational interactions with companion bodies and the low order spin-orbit resonance of the librating planets we 
are considering, the general figure of the CMB must be ellipsoidal with a polar flattening and a tidal bulge pointing 
on average toward the main gravitational partner [15] . 

The fluid dynamics that occurs in a rapidly rotating ellipsoidal cavity has been widely studied in the case of 
constant but different rotation rates of the fluid and the elliptical distortion. This corresponds in geophysical terms 
to a non-synchronized body with a constant spin rate subject to dynamical tides rotating at the constant orbital 
rotation rate f2 or fc. In particular, it has been shown that this elliptically deformed base flow can be destabilized by 
the so-called tidally-driven elliptical instability or TDEI [TBI [T7] . 

Generally speaking, the elliptical instability can be seen as the inherent local instability of elliptical streamlines 
[18H20] . or as the parametric resonance between two free inertial waves (resp. modes) of the rotating unbounded (resp. 
bounded) fluid and an elliptical strain, which is not an inertial wave or mode [5TJ [55] . Such a resonance mechanism, 
confirmed by numerous works in elliptically deformed cylinders [e.g. 23-25 and ellipsoids [26 28 , also operates for 
triadic resonance of three inertial modes, proposed to explain the secondary instability of the elliptical instability 
[29] [30] . This triadic resonance of three inertial modes also applies for the inertial precessional instability in cylinders 
[3TJ[35] and spheroids [531 [33]: there, the base flow forced by precession is itself an inertial mode (e.g. the so-called 
Poincare or tilt-over mode in the spheroid [35 - 38 ), which resonates with two inertial modes. 

It has been shown that selected resonances of the TDEI are sensitive to the ratio of the rotation rates of the 
fluid and the elliptical distortion [28, 39 . In particular, the elliptical instability is known to vanish in the case of 
synchronous rotation fl = Q or i,. But theoretical arguments suggest that oscillations around this synchronous state 
may be sufficient to excite elliptical instability [T71 HOI WT\ . This could be of fundamental importance in planetary 
liquid cores and subsurface oceans of synchronized bodies, where librations are generically present. This letter thus 
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FIG. 1: (a) Mean rotation rates involved in our modelized planetary fluid layer: the elliptical (tidal) deformation rotates at 
the mean orbital rate Q rb, the external solid boundary has a constant tangential velocity imposed by the mean planetary spin 
rate fid and the fluid mean rotation rate in the bulk of the fluid is fio (possibly different from fid in a general case ' L 12l). (b) 
Schematic view of the oscillating triaxial ellipsoidal container. 



aims at validating the existence of a libration driven elliptical instability, hereafter referred to as LDEI. To do so, 
we first extend previous analytical studies [T71 HOI EI] using a local WKB (Wentzel-Kramers-Brillouin) approach that 
allows us to determine a generic formula for the growth rate of LDEI. We then present the numerical and experimental 
validation of the existence of the LDEI, in good agreement with the theoretical results. Finally, implications for planets 
and moons are briefly discussed. 

We consider a homogeneous, non-conductive, incompressible, newtonian fluid enclosed in a librating triaxial ellipsoid 
(fig[l]). In the container frame of reference, the equation of the ellipsoidal boundary is written as x 2 /a 2 +y 2 /b 2 +z 2 jc? = 
1, where (x, y, z) is a cartesian coordinate system with its origin at the center of the ellipsoid, with x along the long 
equatorial axis a, y along the short equatorial axis 6, and z along the rotation axis c. We define the ellipticity 
(3 = (a 2 — b 2 )/(a 2 + b 2 ) and the aspect ratio c* = c/R, where R — J (a 2 + b 2 )/2 is the mean equatorial radius. 
The motion of longitudinal libration of the container is modeled by a time dependency of its axial rotation rate 
r2(t) = Qd + A(j)UJi sin(wjt), where Qd is the mean rotation rate of the container (d for diurnal), A(f> the amplitude of 
libration in radians and uji the angular frequency of libration. In the frame of reference attached to the container, the 
equations of motion, made dimensionless using R as the length scale and Q^ 1 as the time scale, become: 



In ([T]), IT is the reduced pressure, which includes the time- variable centrifugal acceleration. The Ekman number E 
is defined by E = v/(QdR 2 ), where v is the kinematic viscosity, the dimensionless libration frequency / is defined 
as / = uji/Qdi and e is the libration forcing parameter defined by e = A<fif. In the limit of small Ekman numbers, 
the flow can be decomposed into an inviscid component U in the volume and a viscous boundary layer flow u that 
satisfies the no-slip boundary condition on the CMB. Introducing this separation, [3D] proposed the following solution 



du 

~dt 



u x (V x u) + 2 [1 + e sin(/t)] z x u 



VII + E V 2 u - ef cos(/t)(z x r), 



(1) 
(2) 



V-u = 0. 
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to the inviscid equations of motion subject to the non-penetration condition at the CMB: 

U = -Esin(/t) [z x r-/3 V(xy)}. (3) 

It can be shown that eq. Q provides a solution of the inviscid form of ([I]) in the bulk; the inertial forces are 
balanced by the pressure gradient. No net zonal flow can result from the non-linear interactions in the quasi-inviscid 
interior |13j . since pressure gradients alone cannot drive net zonal flows. However, the no-slip boundary condition is 
not entirely fulfilled by this solution: viscous corrections in the Ekman boundary layer must be considered, whose 
non- linear interactions generate a zonal flow in the bulk [TUJ EH [H] , as observed in axisymmetric containers [5HT2"] . 

In addition to these laminar and mostly two-dimensional (2D) motions, Kerswell and Malkus [3D] first suggested 
that turbulent three-dimensional (3D) motions can be excited by an elliptical instability corresponding to a parametric 
resonance involving two free inertial waves of the rotating fluid and strain of the inviscid base flow ([3]) , from which 
energy is extracted. Since the base flow is of azimuthal wavenumber m — 2 and temporal frequency /, this parametric 
resonance occurs only when m a — mj, = ±m = ±2 and A a — At, = ±/, where m a , mb are the azimuthal wavenumbers 
of the two resonating free inertial waves and A a , Ah are their frequencies non-dimensionalized by Qd- In addition to 
these resonance conditions, the two waves must have close radial and azimuthal structures to interact positively, 
corresponding to the so-called principal resonances [32]. This set of rules forms the basis for a global analysis of 
the elliptical instability, which thus requires an exact description of the inertial modes in the considered ellipsoidal 
geometry. Unfortunately, little is known about inertial modes for the finite values of j3 considered in the present study. 
This makes a complete characterization of the elliptical instability by global analysis presently out of reach for our 
simulations and experiments. 

As an alternative, a local stability analysis of the base flow ( |3|) ca n be performed, independently of the geometry 
of the container. Here, we make use of the results presented in |41] for the special case / = I, later complemented 
[T7] . The local analysis is based on the WKB method |20| , which allows an upper bound to be derived for the growth 
rate of the elliptical instability. In this approach, perturbations are assumed to be sufficiently localized so as to be 
advected along flow trajectories. The perturbations are sought in the form of local plane waves characterized by their 
wavevector k(t), with norm k 3> 1, and tilted by an angle £ to the rotation axis. Elliptical instability appears by 
resonance of two identical plane waves, only differing by their direction of propagation. The inviscid growth rate, 
(Jinv, can be determined by solving the Euler equations at the first order in the forcing amplitude fie, yielding 



<Tin» = 16 6 4 /reS ^, (4) 

where f res ^ stands for a resonant forcing frequency [see details in ref. [17]. In the strict WBK limit k ^> 1, the 
dispersion relation between the forcing frequency and the excited plane waves is / res /2 = ±2cos£; hence, all forcing 
frequencies between —4 and +4 should be resonant. However, by accounting for the shape of the container and for 
viscous damping, similarly to the TDEI, resonances are only possible for selected couples of inertial waves, especially 
at the rather large Ekman numbers accessible to numerics [20, 28] . Thus, the system resonates only for selected values 
of the forcing frequency. Introducing a small detuning between the libration frequency / and a given exact resonance 
f res [see method in ref. I2Q[ for the standard case of TDEI] and taking into account the dominant viscous damping in 
the Ekman layer, it can be shown that excitation of the instability takes place around each resonant frequency in a 
band / e [fres — &inv\ fres + &inv\, where the typical growth rate is 
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af nv -(f res -f) 2 -KE^, (5) 



with K € [1; 10] a constant of order 1, which depends on the excited mode. The first term on the left hand side of 
^ defines the range of unstable frequencies around a particular resonance; the second term describes the viscous 
damping of the instability. Besides the strict WBK limit, we expect this equation to be generally valid, once the 
specific values of the resonant frequency f res and of the damping factor K have been determined [see validation for 
the case of TDEI in ref. [2H]- Since both f res and K depend on the selected inertial waves, both values can vary with 
the shape of the container. 

To quantitatively validate the existence of LDEI and the above analytical results, we combine a systematic numerical 
study with selected laboratory experiments. We use the software Comsol Multiphysics® based on the finite elements 
method to perform our simulations. We work in the container frame of reference where the ellipsoidal shape of the 
container is stationary, and we solve the equations of motion ([T]) subject to no-slip boundary conditions, starting from 
a fluid at rest at time t = 0. We refer the reader to [43] for more informations on the numerical method. A first series 
of computations has been performed in 2D to test the realization of the inviscid base flow ([3]). A typical result is 
shown in figure [2] After a transient behavior scaling as a typical viscous time in t ~ E^ 1 (figure^), the theoretical 
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FIG. 2: (color online) (a) Temporal evolution ofU = |u| at the location (x = 0.5, y — 0), as calculated by the 2D version of our 
numerical model for f = 0.5, e = 0.1, (3 — 0.1 and E = 4 x 10 -4 . (b) Zoom (t > 1965,) of figure (a) once the flow is established: 
the continuous gray curve (red online) stands for the numerical results, and the dashed black curve for the theoretical flow pi). 




FIG. 3: (color online) (a) Time evolution of the absolute value of the axial velocity integrated over the whole container, 
non-dimensionalized by the mean value between times t = 100 and t = 200 (i.e. after the spin-up and before the potential 
destabilization) . Numerical simulations are performed for E — 5 x 10 -4 , (3 = 0.44, / = 1.76, £ = 0.92 and c* = 0.95. (b) |u|, 
in a meridional cross section and in the plane z = —0.5. The sequence shows, from left to right, the typical field |u| during the 
exponential growth, at saturation and during the collapse. 



base flow (|3| is indeed established in the bulk (figure |2p) , whereas the corrections necessary to fulfill the tangential 
part of the boundary conditions are restricted to the Ekman layer of depth E~ x l 2 . In the 3D ellipsoid, the Ekman 
pumping associated with the Ekman layer acts to significantly accelerate the establishment of the base flow, and we 
thus expect ^ to be the starting velocity field. 

Series of numerical simulations have then been performed in a triaxial ellipsoid as a function of the libration 
frequency /, the libration amplitude e and the aspect ratio c* , keeping E = 5 x 10~ 4 and (3 — 0.44. An example of 
the temporal evolution of the absolute value of the axial velocity integrated over the volume, W, is presented in figure 
[3^i for the case / = 1.76, e — 0.92 and c* = 0.95. Libration driven elliptical instability is present, characterized by 
intense 3D motions with rich temporal dynamics on typical times much longer than the spin and libration periods. 
In particular, we observe characteristic cycles of growth, saturation, collapses and relaminarization towards the base 
flow, already observed for the classical case of TDEI [3H]. The typical changes in the flow field during one cycle are 
illustrated in figure which shows the norm of the velocity in meridional and equatorial cross sections. 

The growth rate of LDEI can be obtained by fitting the growing parts of time series of the integrated axial velocity 
W with an exponential function. The systematic evolution of this growth rate with / for c* = 1 is shown in figure [4^,, 
in comparison with the analytical formula where f res and K have been determined by adjusting ^ to each of the 
two local maxima of the numerically determined growth rate. We observe two bands of frequency centered around 
f res ~ 1.83 and f res ~ 1.67 where an elliptical instability grows. Good agreement is recovered for all neighboring 
values of /, validating the generic equation 
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FIG. 4: (color online) Systematic numerical determination of the growth rate o of the LDEI for E — 5 x 10 -4 , P = 0.44, 
e — 0.991 (a is set to when the LDEI is not excited), (a) As a function of the libration frequency f for c* = 1. Also shown 
here as a solid curve is the theoretical growth rate where f res and K have been determined by adjusting the analytical 
formula to each of the two local maxima of the numerical growth rate, (b) As a function of the aspect ratio c* for f = 1.8. Note 
that the LDEI is excited in spheroidal geometries (c — b or c = a). 



The systematic evolution of the growth rate as a function of the aspect ratio c* at a hxed resonant frequency 
/ = 1.8 is shown in figure [4]d. Although less dramatically than the frequency detuning, the geometrical factor c* can 
also alter the growth rate of the elliptical instability at a given frequency. We observe that the optimal geometry at 
the resonant frequency / = 1.83 is realized for c* ~ 1, as already shown in [43] for the classical elliptical instability 
(TDEI). In contrast with the TDEI studied in [S3] , we have no theoretical arguments that precludes the excitation 
of the LDEI in spheroidal geometries (c = b or c = a). In figure ^p, we label these two particular cases, showing that 
the growth rates in these geometries are positive. Note that this does not refute the conclusions of [44] . who show 
that libration in longitude cannot produce a direct resonance of an inertial mode. Here, the resonance docs not occur 
directly between a mode and the forcing but between two modes and the forcing in a parametric coupling. The modes 
excited in the LDEI are not at the frequency of the forcing (which would be the case for a direct forcing) , it is their 
frequency difference that is equal to /. 

To further validate the existence of LDEI, we have also performed selected laboratory experiments (fig. [5]). These 
allow us to reach smaller values of the Ekman number, and hence to access more chaotic flows, with the inconvenience 
being the difficulty in acquiring precise quantitative measurements and in changing the shape of the container. Except 
for the container, the laboratory apparatus is the same as in ref. [8] and [TTj . It consists in an oscillating tank filled 
with water and centered on a turntable rotating at a constant angular velocity £7^ = 0.5 Hz. In order to perform 
quantitative measurements using Laser Doppler Velocimetry (LDV), the container consists in one hemispheroid CNC 
(Computer Numerical Control) machined from cast acrylic cylindrical blocks that is polished optically clear, with 
a top flat lid that avoids optical distortions. This geometry, corresponding to a "half-ellipsoidal" container, does 
not allow modes of the elliptical instability that are antisymmetric around the equator. Note also that because 
of manufacturing constraints, the small axis and the rotation axis of the container are equal, with a = 127 mm, 
b = c = 89 mm. The experimental parameters are then j3 — 0.34, c* = 0.812 and E = 2.7 x 10~ 5 ; / and e have been 
changed systematically to explore the ranges / € [0.5 — 2] with e = 1, and e E [0 — 1.6] with / = 1, respectively. 
Similar to the numerical experiments, we observe a resonance band, here / G [1.43 — 1.66], that is characterized by 
strong, intermittent, space-filling turbulence. These cases are marked by periods dominated by strong, small-scale, 
shear structures, followed by relaminarization period. These flows are not due to shear instabilities since the system 
does not exhibit any turbulence at higher frequency, for which the Rossby number is larger. Since no direct resonance 
can occur, if a shear instability were developing for a critical Rossby number it should remain unstable at larger 
forcing, regardless of the frequency. 

An example of the measured azimuthal velocity is shown in figure [6] for / = 1.46. As mentioned earlier, the observed 
chaotic behavior of the flow as well as the fact that this behavior only appears for specific frequencies are characteristic 
of the elliptical instability. Assuming the existence of a resonant peak of LDEI at / = 1.46, the WKB approach ^ 
yields a typical growth time 1 / (crO^) of the instability ranging between 21s for a lower bound of the damping factor 
K = 1 and 37s for an upper bound of the damping factor K = 10. As shown in figure [6j these values are in good 
agreement with measurements during each growing phase of the azimuthal flow, further validating the interpretation 
of this chaotic behavior in terms of LDEI. Finally, note that even if no resonant cases have been obtained in this work 
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FIG. 5: Schematic view of the laboratory experiment and of the ultraLDV system \1 l'j . oriented to perform measurements of 
the azimuthal velocity. No LDV measurements of the azimuthal velocity can be readily performed in the case of an ellipsoidal 
container with a non-axisymmetric equatorial cross-section. Indeed, except if the LDV is moving with the tank, no relative 
position of the LDV with the container can be found such that the two laser beams remain coplanar over a libration cycle. A 
co-librating configuration is not achievable in the present experimental set-up due to the large centrifugal acceleration involved. 
Therefore, LDV measurements are performed with a "half-ellipsoidal" container, the northern half of the ellipsoid being replaced 
with a flat plate of acrylic. 



with / < 1 [the relevant planetary range for / according to ref. H] , the WKB analysis predicts unambiguously that 
the LDEI can be excited for an arbitrary |/| < 4 provided that E is sufficiently small. Experiments at lower E should 
confirm it in the future. But the first experimental results presented here, in addition to the numerical simulations 
and in agreement with the theoretical analysis, already illustrate the generic feature of the libration driven elliptical 
instability, which appears for different geometries and various ranges of parameters. 

In conclusion, studies of flows driven by longitudinal libration made using a spherical CMB approximation suggest 
that purely viscous coupling does not lead to significant planetary energy dissipation, angular momentum transfer nor 
magnetic field induction [5]. These conclusions should be re-addressed in accounting for the specific triaxial shape of 
the considered planets, since space-filling turbulence is observed in the present numerical and laboratory experiments 
in which LDEI is excited. A complete understanding of the elliptical instability excited by libration in celestial bodies 
requires the characterization of the inertial modes, their frequency and their viscous decay factor in the appropriate 
geometry and for the low values of the Ekman number relevant to planetary applications. Nevertheless, the relevance 
of a LDEI mechanism at planetary settings can be ascertained using the theoretical WKB approach presented here, 
by estimating a lower bound of the equatorial ellipticity leading to a positive growth rate a. Using typical values 
for forced longitudinal libration [5] of / = 1, E ~ 10~ 14 , e ~ 10 -4 , a decay factor K = 1 and assuming a perfect 
resonance, equation ^ yields to a minimum equatorial ellipticity /3 ~ 10~ 3 for excitation of LDEI. This first order 
estimate is qualitatively comparable with the values expected for Mercury, Europa or Io. Thus, the space-filling 
turbulence resulting from LDEI should exist within the fluid interiors of librating terrestrial bodies. 
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FIG. 6: (color online) Time evolution of the norm of the azimuthal velocity averaged over 10 oscillations for A(j> = 0.7 rad 
(e ~ 1) with P = 0.34 and f = 1.46 (continuous upper curve, red online); /3 — 0.06 and f = 1.40 (continuous upper curve, 
blue online). The measurements are performed at a cylindrical radius Si = 48 mm along the short axis of the mean equatorial 
ellipse, 1 cm below the top flat surface. We perform a sliding window averaging over 10 oscillations with an overlap of 90%. 
In addition we represent the WKB exponential growth for to two extreme values of the damping factor, K = 1 (dotted black) 
and K = 10 (dashed black). The letters L and T stand for Laminar and Turbulent. The periods of turbulence, as observed in 
direct visualizations, are qualitatively represented by the bands. 
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